Présentation du document

NB : Ce script est un document de travail.

Il cherche à objectiver, par l’utilisation de méthodes d’analyse spatiale, l’« effet de grappe » repéré par la littérature (Bideau 2019)

Il est mis à disposition dans une logique de science ouverte.

Ce travail s’inscrit dans le cadre d’une étude plus générale sur les communes nouvelles (et plus particulièrement d’un addendum au chapitre IV de la thèse) :

https://cv.hal.science/gabriel-bideau

Licence CC-BY-NC-SA.

Il est possible d’accéder au code de ce Markdown ici : https://gbideau.github.io/CN_effet_de_grappe/effet_de_grappe.Rmd

Les données utilisées pour jouer le code sont regroupées ici : https://gbideau.github.io/CN_data/

Ne pas hésiter à contacter l’auteur (gabriel.bideau@gmail.com) pour toute question.

Packages nécessaires

# Librairies utilisées
library(sf)
library(cartography)
library(mapsf)
library(corrplot)
library(cowplot)
library(MTA)
library(readxl)
library(ggplot2)
library(FactoMineR) 
library(factoextra)
library(cluster)
library(reshape)
library(reshape2)
library(flows)
# NB : Pour le package flows, la version la plus récente est disponible ici :
# remotes::install_github("rCarto/flows") # ou # install.packages("mapsf")
# Pour obtenir une version plus ancienne (celle utilisée ici) : https://cran.r-project.org/src/contrib/Archive/flows/
# install.packages("packages/flows_1.1.1.tar.gz", repos=NULL, type="source")
library(sp)
library(knitr)
library(condformat) # https://cran.r-project.org/web/packages/condformat/vignettes/introduction.html
library(units)
library(stringr)
# library(dplyr)
library(questionr)
library(spdep) # Pour les matrices de contiguïté
library(rgeoda) # Pour les matrices de contiguïté

# Liste pour installer les packages si besoin :
# sf cartography mapsf readxl foreign dplyr flextable knitr stringr units condformat forcats ggplot2 rstatix questionr corrplot gtsummary broom GGally effects forestmodel ggeffects labelled cowplot spdep rgeoda

Matrice de contiguïté et autocorrélation spatiale

Import des données

load("data/refdata.Rdata")

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE)

Il s’agit ici de creuser l’« effet de grappe », ce qu’on appelle parfois l’« effet consultant ».

Matrice de voisinage (package spdep)

À l’aide de Bellefon, Loonis, and Le Gleut (2018).

library(spdep)

test_geom <- merge(geom2011, df2011[, c("CODGEO", "CODE_DEPT", "COM_NOUV", "P09_POP", "LIBGEO")], by = "CODGEO")
test_geom <- subset (test_geom, CODE_DEPT == "49" )

test_queen <- poly2nb(test_geom, queen = TRUE)
test_rook <- poly2nb(test_geom, queen = FALSE)

# On change de type d'objet pour faciliter les représentations
test_geom <- as(test_geom, "Spatial")

# Représentation graphique des deux manières de calculer les voisinages.
plot(test_geom, border="lightgray")
plot(test_queen, coordinates(test_geom),add=TRUE,col="red")
plot(test_rook, coordinates(test_geom),add=TRUE,col="blue")

# Réalisation d'une liste de poids
test_liste <- nb2listw(test_queen, zero.policy = TRUE)
# NB : zero.policy :default NULL, use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors https://r-spatial.github.io/spdep/reference/nb2listw.html

# Pour avoir une représentation graphique des communes nouvelles
test_geom2 <- merge(geom2011, df2011[, c("CODGEO", "CODE_DEPT", "COM_NOUV", "P09_POP", "LIBGEO")], by = "CODGEO")
test_geom2 <- subset (test_geom2, CODE_DEPT == "49" )

plot(subset(test_geom2[, c("CODE_DEPT", "COM_NOUV")], CODE_DEPT == "49"))

# Pour extraire les entités sans liens :
test_liste$neighbours

## Neighbour list object:
## Number of regions: 363 
## Number of nonzero links: 2012 
## Percentage nonzero weights: 1.526915 
## Average number of links: 5.5427

isol <- test_geom2[c(5896, 7662, 10958, 10959, 10960, 11024, 21454, 21470, 21471, 21472, 21473, 33974),]


rm(isol, test_queen, test_rook)

À l’échelle du département du Maine-et-Loire, des différences notables en fonction du voisinage choisi.

Officiellement, le contact par un point autorise à fusionner en commune nouvelle donc choix du voisinage “queen”.

Cf. la réponse du ministère de l’intérieur à une question écrite au Sénat sur la question de la continuité territoriale en 2009 : https://www.senat.fr/questions/base/2009/qSEQ090609011.html.

Calcul de l’indice de Moran (package spdep)

Cf. Bellefon and Salima (2018) p. 56 sq. pour l’explication et 62 pour le code.

test_geom$COM_NOUV2 <- dplyr::if_else(test_geom$COM_NOUV == "OUI", 1, 0)

moran.plot(test_geom$COM_NOUV2, test_liste, labels=FALSE , xlab="variable",ylab="moyenne des voisins")

moran.plot(test_geom$P09_POP, test_liste, labels=FALSE , xlab="variable",ylab="moyenne des voisins")

Statistiques join count (package spdep)

Cf. Bellefon and Salima (2018) p. 62 sq. pour l’explication et 64 pour le code.

L’analyse des statistiques des join count observe, pour une variable binaire, s’il y a une autocorrélation spatiale de la présence d’individus partageant ou non une caractéristique. “on considère une variable binaire qui représente deux couleurs, Blanc (B) et Noir (N) de sorte qu’une liaison puisse être qualifiée de Blanc-Blanc, Noir-Noir ou Blanc-Noir. On observe :

— une autocorrélation spatiale positive si le nombre de liaisons Blanc-Noir est significativement inférieur à ce que l’on aurait obtenu à partir d’une répartition spatiale aléatoire ;

— une autocorrélation spatiale négative si le nombre de liaisons Blanc-Noir est significativement supérieur à ce que l’on aurait obtenu à partir d’une répartition spatiale aléatoire ;

— aucune autocorrélation spatiale si le nombre de liaisons Blanc-Noir est approximativement identique à ce que l’on aurait obtenu à partir d’une répartition spatiale aléatoire.” (p. 62-63, souligné dans le texte)

test_geom$COM_NOUV3 <- as.factor(test_geom$COM_NOUV)
class(test_geom$COM_NOUV3)

## [1] "factor"

joincount.test(test_geom$COM_NOUV3, test_liste, zero.policy=TRUE, alternative="greater",
 sampling="nonfree", spChk=NULL, adjust.n=TRUE)

## 
##  Join count test under nonfree sampling
## 
## data:  test_geom$COM_NOUV3 
## weights: test_liste 
## 
## Std. deviate for NON = 12.114, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             44.663005             26.113260              2.344626 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  test_geom$COM_NOUV3 
## weights: test_liste 
## 
## Std. deviate for OUI = 11.408, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             89.519769             69.613260              3.045139

Pour les communes nouvelles, on observe une autocorrélation positive très nette et significative. L’autocorrélation est d’ailleurs également significtive pour les communes non fusionnantes, même si la déviation est moins nette.

LISA : Indice de Moran (package spdep)

Mesures locales d’association spatiale : LISA (Local Indicators of Spatial Association).

Cf. Bellefon and Salima (2018) p. 65 sq. pour l’explication et 69 pour le code.

“Anselin (ANSELIN 1995) développe les notions introduites par Getis et Ord en définissant des indicateurs d’autocorrélation spatiale locale. Ceux-ci doivent permettre de mesurer l’intensité et la significativité de la dépendance locale entre la valeur d’une variable dans une unité spatiale et les valeurs de cette même variable dans les unités spatiales environnantes.” (p. 65)

Fonction ‘localmoran’ du package ‘spdep’. Ii : local moran statistic ; E.Ii : expectation of local moran statistic ; Var.Ii : variance of local moran statistic ; Z.Ii : standard deviate of local moran statistic ; Pr() : p-value of local moran statistic (https://www.rdocumentation.org/packages/spdep/versions/1.1-8/topics/localmoran).

Pour analyser correctement la représentativité des résultats, plusieurs méthodes proposent d’ajuster la p-value :

“On considère que cette méthode [celle de Bonferroni] ne donne de bons résultats que lorsque le nombre de tests réalisés est petit. […] La méthode d’ajustement de Holm conduit à un plus grand nombre de clusters significatifs que la méthode de Bonferroni. Elle lui est donc le plus souvent préférée. Cependant, cette méthode se concentre aussi sur la détection de l’existence d’au moins un cluster dans toute la zone. […] La méthode du False Discovery Rate (FDR) a été introduite par BENJAMINI et al. 1995. Avec cette méthode, le risque de juger - à tort - un cluster comme significatif est plus élevé, mais inversement le risque de juger - à tort - un cluster comme non significatif est plus faible.” (p. 68) “La méthode de Holm diminue en effet le risque de conclure à tort à l’existence d’une autocorrélation spatiale locale. En revanche, cette méthode augmente le risque de passer à côté d’un cluster local.” (p.69)

results_LISA<- as.data.frame(localmoran(test_geom$COM_NOUV2,test_liste,zero.policy=TRUE))

#Calcul des p-values ajustées
results_LISA$pvalue_ajuste_bonferroni<- p.adjust(results_LISA$`Pr(z != E(Ii))`,method="bonferroni")
results_LISA$pvalue_ajuste_holm<- p.adjust(results_LISA$`Pr(z != E(Ii))`,method="holm")
results_LISA$pvalue_ajuste_fdr<- p.adjust(results_LISA$`Pr(z != E(Ii))`,method="fdr")

test_geom3 <- cbind(test_geom2, results_LISA)


# Quelles communes ont une p-value élevée selon les trois méthodes ?
par(mar = c(0,0,0,0), mfrow = c(2,2))

# plot(st_geometry(dep), col = NA)
choroLayer(x = test_geom3 , var = "pvalue_ajuste_holm",
           method = "quantile", nclass = 6,
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = carto.pal(pal1 = "blue.pal", n1 = 6),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Pvalue ajustée : méthode de Holm")


choroLayer(x = test_geom3 , var = "pvalue_ajuste_bonferroni",
           method = "quantile", nclass = 6,
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = carto.pal(pal1 = "blue.pal", n1 = 6),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Pvalue ajustée : méthode de Bonferroni")

choroLayer(x = test_geom3 , var = "pvalue_ajuste_fdr",
           method = "quantile", nclass = 6,
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = carto.pal(pal1 = "blue.pal", n1 = 6),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Pvalue ajustée : méthode FDR")


typoLayer(x = test_geom3, var = "COM_NOUV",  
          col=c("red", "blue"), border = NA)

layoutLayer(title = "Titre", theme = "red.pal",
            author = "G. Bideau, 2022.",
            sources = "")

# On choisit la méthode FDR et on ne conserve que les valeurs qui nous paraissent significatives au vu de la p-value ainsi ajustée

test_geom3$Ii_retenu <- ifelse(test_geom3$pvalue_ajuste_fdr <= 0.1, yes = test_geom3$Ii, no = NA)

par(mar = c(0,0,0,0), mfrow = c(1,2))

# On cartographie le résultat de l'indice de Moran

choroLayer(x = test_geom3 , var = "Ii_retenu",
           method = "quantile", nclass = 6,
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = carto.pal(pal1 = "blue.pal", n1 = 6),
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Indices de Moran significatifs",
           legend.nodata = "p-value non significatif (>0,1 méthode FDR)")
typoLayer(x = test_geom3, var = "COM_NOUV",  
          col=c("red", "blue"), border = NA, legend.pos = "topleft")

par(mar = c(0,0,0,0), mfrow = c(1,1))

Matrice de voisinage (package rgeoda)

Présentation du package ‘rgeoda’ : https://geodacenter.github.io/rgeoda/articles/rgeoda_tutorial.html

Sur les matrices de contiguïté : https://geodacenter.github.io/workbook/4a_contig_weights/lab4a.html

Sur l’auto-corrélation : http://geodacenter.github.io/workbook/6a_local_auto/lab6a.html

library(rgeoda)


test_geom <- merge(geom2011, df2011[, c("CODGEO", "CODE_DEPT", "COM_NOUV", "P09_POP")], by = "CODGEO")
# test_geom <- subset (test_geom, CODE_DEPT == "49" )
test_geomCN <- merge(geom_new, df_new, by = "CODGEO_new")
test_geomCN <- subset(test_geomCN, COM_NOUV == "OUI")

queen_w <- queen_weights(test_geom, order=1, include_lower_order = FALSE, precision_threshold = 0)
summary(queen_w)

##                      name                value
## 1 number of observations:                36208
## 2          is symmetric:                  TRUE
## 3               sparsity: 0.000164767983149941
## 4        # min neighbors:                    0
## 5        # max neighbors:                   29
## 6       # mean neighbors:     5.96591913389306
## 7     # median neighbors:                    6
## 8           has isolates:                 TRUE

# Pour savoir si certaines entités sont séparées des autres
has_isolates(queen_w)

## [1] TRUE

LISA (package rgeoda)

NB : Sous-section qui peut poser problème lors d’un knit si n’est pas executée en même temps que les précédentes (cela ne fonctionne pas bien si les sous-sections précédentes sont en cache).

test_geom$COM_NOUV2 <- dplyr::if_else(test_geom$COM_NOUV == "OUI", 1, 0)
variable <- "COM_NOUV2"
df_variable <- test_geom[variable]
lisa <- local_moran(queen_w, df_variable)

# To get the values of the local Moran:
lms <- lisa_values(gda_lisa = lisa)

# To get the pseudo-p values of significance of local Moran computation:
pvals <- lisa_pvalues(lisa)

# To get the cluster indicators of local Moran computation:
cats <- lisa_clusters(lisa, cutoff = 0.05)

# Labels
lbls <- lisa_labels(lisa)


test_geom$LISA <- lms
test_geom$pvals <- pvals
test_geom$cats <- cats

# On rend les catégories plus lisibles en remplaçant le numéro par l'intitulé
test_geom$cats <- as.factor(test_geom$cats)
num_levels <- as.numeric(levels(test_geom$cats))
levels(test_geom$cats) <- lbls[num_levels + 1]
print(levels(test_geom$cats))

## [1] "Not significant" "High-High"       "Low-High"        "High-Low"       
## [5] "Isolated"

test_geom$LISA_retenu <- ifelse(test_geom$pvals <= 0.1, yes = test_geom$LISA, no = NA)

par(mar = c(0,0,0,0), mfrow = c(1,2))

couleurs <- c("#f7f7f7", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33")
# Cf. https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6 pour composer les palettes
typoLayer(x = test_geom , var = "cats",
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = couleurs[1:length(levels(test_geom$cats))],
           border = NA,
           legend.pos = "topleft",
           legend.values.order = levels(test_geom$cats),
           legend.title.txt = "Indices de Moran significatifs\nsur la variable de « commune nouvelle »\n(pvalue < 0,05)")


typoLayer(x = test_geom, var = "COM_NOUV",
          legend.values.order = c("OUI", "NON"),
          col=c("red","blue"),
          border = NA,
          legend.pos = "topleft",
          legend.title.txt = "Communes fusionnantes")

# plot(test_geomCN$geometry, col = NA, border = "black", lwd = 1, add = TRUE)

par(mar = c(0,0,0,0), mfrow = c(1,1))

table(test_geom$cats, test_geom$COM_NOUV)

##                  
##                     NON   OUI
##   Not significant 31245   992
##   High-High           0  1576
##   Low-High         2380     0
##   High-Low            0     3
##   Isolated           12     0

rm(variable, df_variable, lms, pvals, cats, lbls, num_levels, couleurs)

LISA (package rgeoda) autres variables

L’objectif de cette sous-section est d’observer si on observe une différence entre communes fusionnantes et autres communes, du point de vue de l’auto-corrélation spatiale de certains indicateurs.

Globalement, on observe bien des positionnements parfois un peu différents des communes fusionnantes, mais cela est sans doute davantage lié à des effets de contexte (les communes fusionnantes sont sur-représentés dans certains types d’endroits, sous-représentés dans d’autres) que réellement lié à ces clusters.

test_geom$COM_NOUV2 <- dplyr::if_else(test_geom$COM_NOUV == "OUI", 1, 0)
variable <- "COM_NOUV2"
colnames(df2011)

##  [1] "CODGEO"                 "LIBGEO"                 "CODE_DEPT"             
##  [4] "CATAEU2010"             "REG"                    "ARR"                   
##  [7] "CV"                     "UU2010"                 "AU2010"                
## [10] "ZE2010"                 "EPCI"                   "P09_ACT1564"           
## [13] "P09_CHOM1564"           "P09_ETUD1564"           "P09_RETR1564"          
## [16] "C09_ACT1564_Agr"        "C09_ACT1564_ArtCom"     "C09_ACT1564_Cadr"      
## [19] "C09_ACT1564_ProfInt"    "C09_ACT1564_Empl"       "C09_ACT1564_Ouvr"      
## [22] "P09_EMPLT"              "C09_EMPLT_AGRI"         "C09_EMPLT_INDUS"       
## [25] "C09_EMPLT_CONST"        "C09_EMPLT_CTS"          "C09_EMPLT_APESAS"      
## [28] "P09_ACTOCC"             "P09_POP"                "P09_POP0014"           
## [31] "P09_POP1529"            "P09_POP3044"            "P09_POP4559"           
## [34] "P09_POP6074"            "P09_POP75P"             "C09_ACTOCC_IN"         
## [37] "C09_ACTOCC_OUT"         "C09_ACTOCC"             "P11_POT_FIN"           
## [40] "P11_DGF"                "P11_FoyFisc"            "P11_Rev_Fisc"          
## [43] "P11_IMP_NET"            "P11_FoyFisc_Imp"        "superficie"            
## [46] "CIF_2012"               "CIF_2015"               "ZAU_POL"               
## [49] "ZAU_RUR"                "ZAU_MAR_SP"             "ZAU_MAR"               
## [52] "ZAU_PERI"               "ZAU_AU"                 "FUSION"                
## [55] "ChefLieu"               "ComDLG"                 "FusDate"               
## [58] "FusPhas"                "COM_NOUV"               "CODGEO_new"            
## [61] "LIBGEO_new"             "P09_CHOM1564_RT"        "P09_ETUD1564_RT"       
## [64] "P09_RETR1564_RT"        "C09_ACT1564_Agr_RT"     "C09_ACT1564_ArtCom_RT" 
## [67] "C09_ACT1564_Cadr_RT"    "C09_ACT1564_ProfInt_RT" "C09_ACT1564_Empl_RT"   
## [70] "C09_ACT1564_Ouvr_RT"    "C09_EMPLT_AGRI_RT"      "C09_EMPLT_INDUS_RT"    
## [73] "C09_EMPLT_CONST_RT"     "C09_EMPLT_CTS_RT"       "C09_EMPLT_APESAS_RT"   
## [76] "P09_POP0014Y_RT"        "P09_POP1529Y_RT"        "P09_POP3044Y_RT"       
## [79] "P09_POP4559Y_RT"        "P09_POP6074Y_RT"        "P09_POP75PY_RT"        
## [82] "C09_ACTOCC_IN_RT"       "C09_ACTOCC_OUT_RT"      "C09_EMP_CONC_RT"       
## [85] "P11_POT_FIN_RT"         "P11_DGF_RT"             "P11_Rev_Fisc_RT"       
## [88] "P11_IMP_NET_RT"         "P11_FoyFisc_Imp_RT"

variable <- "P11_POT_FIN"
variable <- "P09_POP1529Y_RT"
variables_a_tester <- c("P09_CHOM1564", "C09_EMPLT_AGRI","P09_POP0014", "P09_POP1529", "P09_POP6074", "C09_ACTOCC_OUT", "superficie", "P11_POT_FIN")

for (variable in variables_a_tester) {

test_geom2 <- merge(test_geom, df2011[, c("CODGEO", variable)], by = "CODGEO")
colnames(test_geom2)
df_variable <- test_geom2[variable]
lisa <- local_moran(queen_w, df_variable)

# To get the values of the local Moran:
lms <- lisa_values(gda_lisa = lisa)

# To get the pseudo-p values of significance of local Moran computation:
pvals <- lisa_pvalues(lisa)

# To get the cluster indicators of local Moran computation:
cats <- lisa_clusters(lisa, cutoff = 0.05)

# Labels
lbls <- lisa_labels(lisa)


test_geom2$LISA <- lms
test_geom2$pvals <- pvals
test_geom2$cats <- cats

# On rend les catégories plus lisibles en remplaçant le numéro par l'intitulé
test_geom2$cats <- as.factor(test_geom2$cats)
num_levels <- as.numeric(levels(test_geom2$cats))
levels(test_geom2$cats) <- lbls[num_levels + 1]
print(levels(test_geom2$cats))

test_geom2$LISA_retenu <- ifelse(test_geom2$pvals <= 0.1, yes = test_geom2$LISA, no = NA)

par(mar = c(0,0,0,0), mfrow = c(1,2))

couleurs <- c("#f7f7f7", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33")
# Cf. https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6 pour composer les palettes
typoLayer(x = test_geom2 , var = "cats",
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = couleurs[1:length(levels(test_geom2$cats))],
           border = NA,
           legend.pos = "topleft",
           legend.values.order = levels(test_geom2$cats),
           legend.title.txt = paste0("Indices de Moran significatifs\nsur la variable « ", variable, " »\n(pvalue < 0,05)"))


typoLayer(x = test_geom2, var = "COM_NOUV",
          legend.values.order = c("OUI", "NON"),
          col=c("red","blue"),
          border = NA,
          legend.pos = "topleft",
          legend.title.txt = "Communes fusionnantes")

# plot(test_geom2CN$geometry, col = NA, border = "black", lwd = 1, add = TRUE)

par(mar = c(0,0,0,0), mfrow = c(1,1))



# On regarde la différence entre les communes fusionnantes et les autres
tabcont <- table(test_geom2$cats, test_geom2$COM_NOUV)

print(tabcont) # En valeur absolue
print(round(100*prop.table(tabcont,margin=1),1)) # Pourcentages, le total se fait par lignes
# round(100*prop.table(tabcont,margin=),1) # Pourcentages, le total se fait sur l'ensemble de la population
print(round(100*prop.table(tabcont,margin=2),1)) # Pourcentages, le total se fait par colonnes

test<-chisq.test(tabcont)
# test$observed
# round(test$expected,1)
# round(test$residuals,2)
print(test)



  
}

## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 25521  2148
##   High-High        1553    27
##   Low-Low          5498   352
##   Low-High          937    30
##   High-Low          116    14
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.2   7.8
##   High-High        98.3   1.7
##   Low-Low          94.0   6.0
##   Low-High         96.9   3.1
##   High-Low         89.2  10.8
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 75.9 83.5
##   High-High        4.6  1.1
##   Low-Low         16.3 13.7
##   Low-High         2.8  1.2
##   High-Low         0.3  0.5
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 125.46, df = 5, p-value < 2.2e-16
## 
## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 24929  1921
##   High-High        2724   294
##   Low-Low          4425   235
##   Low-High         1188    97
##   High-Low          359    24
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.8   7.2
##   High-High        90.3   9.7
##   Low-Low          95.0   5.0
##   Low-High         92.5   7.5
##   High-Low         93.7   6.3
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 74.1 74.7
##   High-High        8.1 11.4
##   Low-Low         13.2  9.1
##   Low-High         3.5  3.8
##   High-Low         1.1  0.9
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 63.652, df = 5, p-value = 2.133e-12
## 
## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 24406  2117
##   High-High        1902    38
##   Low-Low          6576   374
##   Low-High          603    24
##   High-Low          138    18
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.0   8.0
##   High-High        98.0   2.0
##   Low-Low          94.6   5.4
##   Low-High         96.2   3.8
##   High-Low         88.5  11.5
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 72.6 82.3
##   High-High        5.7  1.5
##   Low-Low         19.5 14.5
##   Low-High         1.8  0.9
##   High-Low         0.4  0.7
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 155.87, df = 5, p-value < 2.2e-16
## 
## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 24796  2150
##   High-High        1661    34
##   Low-Low          6315   355
##   Low-High          741    21
##   High-Low          112    11
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.0   8.0
##   High-High        98.0   2.0
##   Low-Low          94.7   5.3
##   Low-High         97.2   2.8
##   High-Low         91.1   8.9
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 73.7 83.6
##   High-High        4.9  1.3
##   Low-Low         18.8 13.8
##   Low-High         2.2  0.8
##   High-Low         0.3  0.4
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 153.54, df = 5, p-value < 2.2e-16
## 
## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 25292  2180
##   High-High        2005    42
##   Low-Low          5495   302
##   Low-High          676    25
##   High-Low          157    22
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.1   7.9
##   High-High        97.9   2.1
##   Low-Low          94.8   5.2
##   Low-High         96.4   3.6
##   High-Low         87.7  12.3
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 75.2 84.8
##   High-High        6.0  1.6
##   Low-Low         16.3 11.7
##   Low-High         2.0  1.0
##   High-Low         0.5  0.9
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 161.05, df = 5, p-value < 2.2e-16
## 
## [1] "High-High" "Isolated"

##            
##               NON   OUI
##   High-High 33625  2571
##   Isolated     12     0
##            
##               NON   OUI
##   High-High  92.9   7.1
##   Isolated  100.0   0.0
##            
##             NON OUI
##   High-High 100 100
##   Isolated    0   0
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabcont
## X-squared = 0.15665, df = 1, p-value = 0.6923
## 
## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 23706  1902
##   High-High        3974   279
##   Low-Low          4376   264
##   Low-High         1228   104
##   High-Low          341    22
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.6   7.4
##   High-High        93.4   6.6
##   Low-Low          94.3   5.7
##   Low-High         92.2   7.8
##   High-Low         93.9   6.1
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 70.5 74.0
##   High-High       11.8 10.9
##   Low-Low         13.0 10.3
##   Low-High         3.7  4.0
##   High-Low         1.0  0.9
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 22.554, df = 5, p-value = 0.0004106
## 
## [1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
## [5] "High-Low"        "Isolated"

##                  
##                     NON   OUI
##   Not significant 25369  2177
##   High-High        1665    32
##   Low-Low          5741   323
##   Low-High          749    27
##   High-Low          101    12
##   Isolated           12     0
##                  
##                     NON   OUI
##   Not significant  92.1   7.9
##   High-High        98.1   1.9
##   Low-Low          94.7   5.3
##   Low-High         96.5   3.5
##   High-Low         89.4  10.6
##   Isolated        100.0   0.0
##                  
##                    NON  OUI
##   Not significant 75.4 84.7
##   High-High        4.9  1.2
##   Low-Low         17.1 12.6
##   Low-High         2.2  1.1
##   High-Low         0.3  0.5
##   Isolated         0.0  0.0
## 
##  Pearson's Chi-squared test
## 
## data:  tabcont
## X-squared = 144.26, df = 5, p-value < 2.2e-16

rm(variable, df_variable, lms, pvals, cats, lbls, num_levels, couleurs)

Cartographies

Sans lien direct avec l’effet de grappe, éléments de cartographie, qui sont développés dans le chapitre X de la thèse (utilisation du même Markdown pour limiter la création de multiples dépôts).

Cf. l’idée donnée ici : https://www.bnsp.insee.fr/ark:/12148/bc6p08tp83z/f1.pdf#page=2

Préparation des données

geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE) 
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) 
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)  
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE)

# Métadonnées
# Liste toutes les variables disponiles
variables_dispo <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ind_target"))
# Liste les variables marquées dans le fichier meta_budgets.xlsx comme nous intéressant
target <- subset(variables_dispo, variable_selec == "X")

Les données socio-économiques qui décrivent les communes en géographies 2011 et 2021 sont ici importées. On commence par extraire les communes ayant participé à la création d’une commune nouvelle, appelées ici communes fusionnantes (datafus2011), les communes nouvelles, avec les géométries au 1er janvier 2021 et caractérisées par les données à la géométrie 2011 agrégées (dataCN_new), ainsi que les communes, à la géométrie 2011, qui n’ont pas participé à la création d’une commune nouvelle (dataNfus2011)

load("data/refdata.Rdata")
datafus2011 <- subset(df2011, COM_NOUV == "OUI")
dataCN_new <- subset(df_new, COM_NOUV == "OUI")
dataNfus2011 <- subset(df2011, COM_NOUV == "NON")

Dans un certain nombre de cas, il sera utile d’avoir, dans un même objet, les données et les géométries. Les données sont ici jointes aux couches géographiques d’intérêt.

geom2011 <- merge(geom2011, df2011, by = "CODGEO")
geom_new <- merge(geom_new, df_new, by = "CODGEO_new")
geomCN_new <- merge(geomCN_new, dataCN_new, by = "CODGEO_new")
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")

En fonction de la taille démographique

# par(mfrow = c(1,1), mar=c(0,0,1.5,0))

# En fonction de la population de la commune fusionnante

choroLayer(x = geomfus2011 , var = "P09_POP",
           method = "quantile", nclass = 4,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 2, n2 = 2),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Nombre d'habitants (2009,\nregroupement par quartiles)")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)


layoutLayer(title = " ",# "Communes fusionnantes (2011-2024)\nen fonction de leur nombre d'habitants",
            author = "G. Bideau, 2024",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2024")

# On crée un tableau donnant, pour chaque département, la population moyenne des communes fusionnantes
tableau_dep <- as.data.frame(tapply(geomfus2011$P09_POP, INDEX = geomfus2011$CODE_DEPT, mean))
colnames(tableau_dep) <- c("moyenne_Cfus")
tableau_dep$median_Cfus <- tapply(geomfus2011$P09_POP, INDEX = geomfus2011$CODE_DEPT, median)
tableau_dep$CODE_DEPT <- row.names(tableau_dep)
tableau_dep$Nbr_Cfus <- table(geomfus2011$CODE_DEP)
tableau_dep <- merge(tableau_dep, dep[, c("CODE_DEPT", "LIBELLE")], by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_dep_tt <- as.data.frame(tapply(geom2011$P09_POP, INDEX = geom2011$CODE_DEPT, mean))
colnames(tableau_dep_tt) <- c("moyenne")
tableau_dep_tt$median <- tapply(geom2011$P09_POP, INDEX = geom2011$CODE_DEPT, median)
tableau_dep_tt$CODE_DEPT <- row.names(tableau_dep_tt)
tableau_dep <- merge(tableau_dep, tableau_dep_tt, by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_pr_publi <- subset(tableau_dep, tableau_dep$Nbr_Cfus > 40)
tableau_pr_publi <- tableau_pr_publi[order(-tableau_pr_publi$moyenne_Cfus),]

kable(tableau_pr_publi[, c("LIBELLE", "Nbr_Cfus", "moyenne_Cfus", "moyenne", "median_Cfus", "median")],
      col.names = c("Département", "Nombre de communes fusionnantes", "Population moyenne des communes fusionnantes", "Population communale moyenne", "Population médiane des communes fusionnantes", "Population communale médiane"), digits = 0)

	Département	Nombre de communes fusionnantes	Population moyenne des communes fusionnantes	Population communale moyenne	Population médiane des communes fusionnantes	Population communale médiane
76	Vendée	44	2455	2221	1283.0	1288.0
45	Maine-et-Loire	225	1339	2149	978.0	1030.0
46	Manche	206	1128	828	379.5	381.0
1	Ain	45	924	1405	351.0	756.0
68	Savoie	48	894	1348	443.5	549.0
73	Deux-Sèvres	73	862	1201	466.0	579.0
70	Seine-Maritime	48	849	1678	474.5	494.0
24	Eure	127	720	863	377.0	416.0
25	Eure-et-Loir	54	655	1056	423.5	430.0
12	Calvados	221	588	964	317.0	340.0
22	Doubs	43	587	884	373.0	255.5
14	Charente	64	579	870	356.5	387.0
21	Dordogne	78	571	740	270.5	359.0
57	Orne	150	474	579	236.0	244.0
80	Yonne	43	392	755	307.0	354.0
36	Jura	78	373	480	159.0	212.5
44	Lozère	50	303	417	179.0	181.0

# En fonction de la population de la commune nouvelle
choroLayer(x = geomCN_new , var = "P09_POP",
           method = "quantile", nclass = 4,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 2, n2 = 2),
           border = NA,
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Nombre d'habitants (2009,\nregroupement par quartiles)")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
plot(st_geometry(geomCN_new), add = TRUE, lwd = 0.2)

layoutLayer(title = "",# "Communes nouvelles (2011-2024)\nen fonction de leur nombre d'habitants",
            author = "G. Bideau, 2024",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2024")

# On crée un tableau donnant, pour chaque département, la population moyenne des communes nouvelles
tableau_dep <- as.data.frame(tapply(geomCN_new$P09_POP, INDEX = geomCN_new$CODE_DEPT_new, mean))
colnames(tableau_dep) <- c("moyenne_CN")
tableau_dep$median_CN <- tapply(geomCN_new$P09_POP, INDEX = geomCN_new$CODE_DEPT_new, median)
tableau_dep$CODE_DEPT <- row.names(tableau_dep)
tableau_dep$Nbr_Cfus <- table(geomfus2011$CODE_DEP) # On garde le nombre de communes fusionnantes pour faciliter les comparaisons entre les deux tableaux
tableau_dep <- merge(tableau_dep, dep[, c("CODE_DEPT", "LIBELLE")], by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_dep_tt <- as.data.frame(tapply(geom_new$P09_POP, INDEX = geom_new$CODE_DEPT_new, mean))
colnames(tableau_dep_tt) <- c("moyenne")
tableau_dep_tt$median <- tapply(geom_new$P09_POP, INDEX = geom_new$CODE_DEPT_new, median)
tableau_dep_tt$CODE_DEPT <- row.names(tableau_dep_tt)
tableau_dep <- merge(tableau_dep, tableau_dep_tt, by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)

tableau_pr_publi <- subset(tableau_dep, tableau_dep$Nbr_Cfus > 40)
tableau_pr_publi <- tableau_pr_publi[order(-tableau_pr_publi$moyenne_CN),]

kable(tableau_pr_publi[, c("LIBELLE", "Nbr_Cfus", "moyenne_CN", "moyenne", "median_CN", "median")],
      col.names = c("Département", "Nombre de communes nouvelles", "Population moyenne des communes nouvelles", "Population communale moyenne", "Population médiane des communes nouvelles", "Population communale médiane"), digits = 0)

	Département	Nombre de communes nouvelles	Population moyenne des communes nouvelles	Population communale moyenne	Population médiane des communes nouvelles	Population communale médiane
45	Maine-et-Loire	225	7922	4431	5631.5	1569.5
76	Vendée	44	6279	2457	3529.0	1421.0
46	Manche	206	4658	1120	2508.5	444.0
70	Seine-Maritime	48	4075	1768	3228.5	503.0
12	Calvados	221	3012	1289	1662.0	407.5
68	Savoie	48	2684	1506	2697.5	612.0
73	Deux-Sèvres	73	2621	1431	1442.5	691.0
24	Eure	127	2473	996	1385.0	444.0
57	Orne	150	2368	759	2138.5	274.0
1	Ain	45	2309	1502	898.5	799.5
25	Eure-et-Loir	54	2210	1166	1522.0	472.0
21	Dordogne	78	1857	819	1216.5	383.0
14	Charente	64	1685	971	1203.5	426.0
80	Yonne	43	1533	812	1067.0	381.0
22	Doubs	43	1401	923	1323.5	260.0
36	Jura	78	1039	529	825.5	235.0
44	Lozère	50	890	508	579.0	209.0

En fonction de la taille démographique, rapportée aux autres communes

On regarde la position des communes fusionnantes puis des communes nouvelles vis-à-vis des communes françaises

par(mfrow = c(1,2), mar=c(0,0,1.2,0))

# Création des données concernant les quantiles, mais en s'appuyant sur les données de l'ensemble des communes françaises
geomfus2011$P09_POP_quantiles <- as.factor(cut(geomfus2011$P09_POP,
                                     breaks=c(quantile(geom2011$P09_POP, probs=seq(0, 1, 1/6), na.rm = TRUE))))
geomCN_new$P09_POP_quantiles <- as.factor(cut(geomCN_new$P09_POP,
                                     breaks=c(quantile(geom_new$P09_POP, probs=seq(0, 1, 1/6), na.rm = TRUE))))

# En fonction de la population de la commune fusionnante

typoLayer(x = geomfus2011 , var = "P09_POP_quantiles",
          border = NA,
          legend.values.order = levels(geomfus2011$P09_POP_quantiles),
          col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
          # col = carto.pal(pal1 = "red.pal", n1 = 5),
          legend.pos = "topleft",
          legend.title.txt = "Situation au sein\ndes quantiles des\ncommunes françaises")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)


layoutLayer(title = "Communes fusionnantes (2011-2024)\nen fonction de leur nombre d'habitant",
            author = "G. Bideau, 2024",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2024")


# En fonction de la population de la commune nouvelle
typoLayer(x = geomCN_new , var = "P09_POP_quantiles",
          col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
          border = NA,
          legend.values.order = levels(geomCN_new$P09_POP_quantiles),
          legend.pos = "topleft",
          legend.title.txt = "Situation au sein\ndes quantiles des\ncommunes françaises")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
plot(st_geometry(geomCN_new), add = TRUE, lwd = 0.2)

layoutLayer(title = "Communes nouvelles (2011-2024)\nen fonction de leur nombre d'habitant",
            author = "G. Bideau, 2024",
            tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
            sources = "Sources : INSEE, IGN, 2024")

On voit bien que les communes nouvelles sont bien davantage dans la strate supérieure. Mais il ne faut pas croire non plus que seules de petites communes fusionnent, ce n’est pas le cas.

Bibliographie

Bellefon, Marie-Pierre de, Vincent Loonis, and Ronan Le Gleut. 2018. “Codifier La Structure de Voisinage.” In Manuel d’analyse Spatiale. Théorie Et Mise En Œuvre Pratique Avec R. Insee Méthodes n\(^\circ\)131, 33–50. https://www.insee.fr/fr/statistiques/fichier/3635442/imet131-f-chapitre-2.pdf.

Bellefon, Marie-Pierre de, and Bouayad-Agha Salima. 2018. “Indices d’autocorrélation Spatiale.” In Manuel d’analyse Spatiale. Théorie Et Mise En Œuvre Pratique Avec R. Insee Méthodes n\(^\circ\)131, 53–72. https://www.insee.fr/fr/statistiques/fichier/3635442/imet131-f-chapitre-2.pdf.

Bideau, Gabriel. 2019. “Les communes nouvelles françaises (2010-2019) : Une réforme territoriale silencieuse.” Annales de Géographie 728 (4/2019): 57–85. https://www.cairn.info/revue-annales-de-geographie-2019-4-page-57.htm.

#CN_effet_de_grappe Objectiver l’effet de grappe concernant les communes nouvelles

G. Bideau

29 septembre 2024 19:23:34